Using RcppWavelet

Installation

library(devtools)
devtools::install_github("msummersgill/RcppWavelet")

Load Libraries

suppressMessages({
  library(RcppWavelet)
  library(data.table)
  library(plotly)
})

Generate Data

set.seed(1234)
sampleLength <- 2000L
sampleFrequency <- 10L

Index <- seq.default(from = 1/sampleFrequency, by = 1/sampleFrequency, length.out = sampleLength)

Scale <- c(seq.default(from = 0, to = 1.25, length.out = sampleLength/2),
           seq.default(from = 1.25, to = 0, length.out = sampleLength/2))

Drift   <- seq.default(from = 0.75, to = 1.5, length.out = sampleLength)

DT <- data.table(Time = Index,
                 c1 = sin(1/8*(2*pi)*Index),
                 c2 = sin(1/Drift*(2*pi)*Index*1.25),
                 c3 = sin(1/5*(2*pi)*Index)*(Scale^3),
                 Noise = rnorm(sampleLength,mean = 0,sd = 0.5))

DT[,Signal := c1 + c2 + c3 + Noise]

Input Signal

DT %>% 
  plot_ly() %>% 
  add_lines(x ~Time, y = ~c1 ,name = "c1 - Constant",
            color = I("goldenrod1"), line = list(width = 1.5)) %>% 
  add_lines(x ~Time, y = ~c2 ,name = "c2 - Shifting Frequency",
            color = I("firebrick2"), line = list(width = 1.5)) %>% 
  add_lines(x ~Time, y = ~c3 ,name = "c3 - Modulated Amplitude",
            color = I("dodgerblue3"), line = list(width = 1.5)) %>% 
  add_lines(x ~Time, y = ~Noise ,name = "Noise",
            color = I("gray30"),line = list(width = 1)) %>% 
  add_lines(x ~Time, y = ~Signal ,name = "Composite Signal",
            color = I("black"), line = list(width = 2.5)) %>% 
  layout(hovermode = "compare",
         plot_bgcolor = "rgba(235,235,235,0.7)",
         paper_bgcolor = "rgba(0,0,0,0)",
         legend = list(x = 0.5, y = -0.1,
                       xanchor = "center",yanchor = "top",
                       orientation = "h",
                       bgcolor = "transparent"),
         yaxis = list(title = "Amplitude",
                      gridcolor = "rgba(255,255,255,1)"),
         xaxis = list(title = "",
                      gridcolor = "rgba(255,255,255,1)"))

Construct a Wavelet Representation

Wavelet <- RcppWavelet::analyze(x = DT[,Signal],
                                bands_per_octave = 64,
                                frequency_min = 1/10,
                                frequency_max = 2,
                                samplerate_hz = sampleFrequency,
                                mother_wavelet = "MORLET",
                                morlet_carrier = 30,
                                cores = 6L)

DT[,Recon := RcppWavelet::reconstruct(wavelet = Wavelet,scale = TRUE,method = "quantile")]

cat(Wavelet$configuration)

## Wavelet Filter:
##  Frequency Range: 0.1 2
##  Bands per Octave: 64
##  Optimisation: 0
## Wavelet:
##  Sampling rate: 10
##  Scale: 0.2
##  Equivalent Frequency (Hz): 23.8865
##  Window Size: 1
##  Type: Morlet
##  Omega0 (carrier frequency): 30

Review Scalogram and Reconstruction

DT %>% 
  plot_ly() %>% 
  add_lines(x ~Time, y = ~Signal ,name = "Original Signal",
            color = I("black"), line = list(width = 1)) %>% 
  add_lines(x ~Time, y = ~Recon ,name = "Wavelet Reconstruction",
            color = I("firebrick2")) %>% 
  layout(hovermode = "compare",
         plot_bgcolor = "rgba(235,235,235,0.7)",
         paper_bgcolor = "rgba(0,0,0,0)",
         legend = list(x = 0.5, y = -0.1,
                       xanchor = "center", yanchor = "top",
                       orientation = "h",
                       bgcolor = "transparent"),
         yaxis = list(title = "Amplitude",
                      gridcolor = "rgba(255,255,255,1)"),
         xaxis = list(title = "",
                      gridcolor = "rgba(255,255,255,1)"))-> ReconstructionPlot

plot_ly() %>% 
  add_heatmap(z = t(Mod(Wavelet$scalogram)),x = Index, y = Wavelet$periods,
              name = "Scalogram") %>%
  layout(yaxis = list(title = "Period"),
         xaxis = list(title = "")) %>% 
  hide_colorbar() -> ScalogramPlot

subplot(ScalogramPlot,ReconstructionPlot, nrows = 2, shareX = TRUE,titleY = TRUE)